% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (c) 2014 Code_FEM developers


function ADD_DOMAIN(node,edge,element)
% Import a mesh generated by PDETOOL as an additional domain

global N_DIM N_NODE NODE ELEMENT_ID
global N_ELEMENT N_NODE_ELEMENT ELEMENT
global N_DOMAIN ELEMENT_DOMAIN NODE_DOMAIN

% If this is the first domain, the procedure is a bit different
if N_DOMAIN == 0
    N_DIM = size(node,1);
    N_DOMAIN = 1;
    % Let's add the nodes
    NODE = node;
    % The number of domain elements
    n_element = size(element,2);
    % The number of boundary elements
    n_edge = size(edge,2);
    % For the domain elements
    % Define the number of nodes per element.
    % We assume that all the elements in this domain have the same number of
    % nodes per element
    N_NODE_ELEMENT = ones(1,n_element)*(size(element,1)-1);
    % Define the list of nodes for each element
    ELEMENT = element(1:end-1,:);
    % Define the description of each element
    ELEMENT_DOMAIN = [N_DOMAIN*ones(1,size(element,2));zeros(1,size(element,2));element(end,:);2*ones(1,size(element,2))];
    % For the boundary elements
    % Add the number of nodes per element.
    N_NODE_ELEMENT = [N_NODE_ELEMENT ones(1,n_edge)*(size(edge,1)-5)];
    % Add the list of nodes for each element
    ELEMENT = [ELEMENT [edge(1:end-5,:)+N_NODE;zeros(size(ELEMENT,1)-size(edge,1)+5,size(edge,2))]];
    % Add the description of each element
    ELEMENT_DOMAIN = [ELEMENT_DOMAIN [N_DOMAIN*ones(1,size(edge,2));ones(1,size(edge,2));edge(end-2,:);ones(1,size(edge,2))]];
    % Total number of elements
    N_ELEMENT = N_ELEMENT + n_element + n_edge;
    % Total number of nodes
    N_NODE = size(node,2);
    % The domain to which each node belongs
    NODE_DOMAIN = ones(1,N_NODE);
else
    % In that case we are adding a domain
    N_DOMAIN = N_DOMAIN + 1;
    % Add the new nodes
    NODE = [NODE node];
    % The number of domain elements
    n_element = size(element,2);
    % The number of boundary elements
    n_edge = size(edge,2);
    % The number of nodes
    n_node = size(node,2);
    % For the domain elements
    % Define the number of nodes per element.
    N_NODE_ELEMENT = [N_NODE_ELEMENT ones(1,n_element)*(size(element,1)-1)];
    % Define the list of nodes for each element
    if size(ELEMENT,1)>=(size(element,1)-1)
        ELEMENT = [ELEMENT [element(1:end-1,:)+N_NODE;zeros(size(ELEMENT,1)-size(element,1)+1,size(element,2))]];
    else
        ELEMENT = [[ELEMENT;zeros(size(element,1)-1-size(ELEMENT,1),size(ELEMENT,2))] element(1:end-1,:)+N_NODE];
    end
    % Define the description of each element
    ELEMENT_DOMAIN = [ELEMENT_DOMAIN [N_DOMAIN*ones(1,size(element,2));zeros(1,size(element,2));element(end,:)]];
    % For the boundary elements
    % Add the number of nodes per element.
    N_NODE_ELEMENT = [N_NODE_ELEMENT ones(1,n_edge)*(size(edge,1)-5)];
    % Add the list of nodes for each element
    ELEMENT = [ELEMENT [edge(1:end-5,:)+N_NODE;zeros(size(ELEMENT,1)-size(edge,1)+5,size(edge,2))]];
    % Add the description of each element
    ELEMENT_DOMAIN = [ELEMENT_DOMAIN [N_DOMAIN*ones(1,size(edge,2));ones(1,size(edge,2));edge(end-2,:)]];
    % Total number of elements
    N_ELEMENT = N_ELEMENT+n_element+n_edge;
    % Total number of nodes
    N_NODE = N_NODE+n_node;
    % The domain to which each node belongs
    NODE_DOMAIN = [NODE_DOMAIN N_DOMAIN*ones(1,n_node)];
end

ELEMENT_ID = 1:N_ELEMENT;
